#if !defined(quasinewton_INCLUDED_)
#define quasinewton_INCLUDED_

#include "nr3.h"

//#define HACK_TO_DIAGNOSIS  // prints out to iteration number versus function value
                           // prints out cif for each iteration in an appended file

#ifdef HACK_TO_DIAGNOSIS
  #include <fstream>
  #include <iostream>
  #include <string>
  #include <iomanip>
  #include <sstream>
  using namespace std;
  static int fileCounter = 1;
  static ofstream qnOut;
#endif

template <class T>
void dfpmin(VecDoub_IO &p, const Doub gtol, Int &iter, Doub &fret, T &funcd, string& directory)
{
  #ifdef HACK_TO_DIAGNOSIS
    char buf [4];
    string name = "qn_func_val";
    qnOut.open( (directory + name + itoa(fileCounter,buf,10) + ".txt" ).c_str() ); 
    qnOut.precision(50);

    ofstream outCif( (directory + "QN_downhill_CIF" + itoa(fileCounter,buf,10) + ".cif").c_str() );
  #endif


	const Int ITMAX=20000;
	const Doub EPS=numeric_limits<Doub>::epsilon();
	const Doub TOLX=4*EPS,STPMX=100.0;
	Bool check;
	Doub den,fac,fad,fae,fp,stpmax,sum=0.0,sumdg,sumxi,temp,test;
	Int n=p.size();
	VecDoub dg(n),g(n),hdg(n),pnew(n),xi(n);
	MatDoub hessin(n,n);
	fp=funcd(p);
	funcd.df(p,g);
	for (Int i=0;i<n;i++) {
		for (Int j=0;j<n;j++) hessin[i][j]=0.0;
		hessin[i][i]=1.0;
		xi[i] = -g[i];
		sum += p[i]*p[i];
	}
	stpmax=STPMX*MAX(sqrt(sum),Doub(n));

  // Perhaps insert a while loop here so that it doesn't max out!
  /***************

    bool end = false;
    while(!end)
    {
      ...
      if (test < TOLX)
      {
        end = true;
      }
      ...
      if (test < gtol)
      {
        end = true;
      }
      ...
    }
    return;
  
  ***************/


  #ifdef HACK_TO_DIAGNOSIS
    qnOut << 0 << "   " << fret << endl;

    stringstream readIter0;
    readIter0 << 0 << "_";
    readIter0 << fret;
    funcd.print_sol(p, fret, outCif, string("Iter_")+readIter0.str());
  #endif

	for (Int its=0;its<ITMAX;its++) {
		iter=its;
    // here: p    - starting point
    //       xi   - search direction
    //       pnew - new point
		lnsrch(p,fp,g,xi,pnew,fret,stpmax,check,funcd);
		fp=fret;
		for (Int i=0;i<n;i++) {
			xi[i]=pnew[i]-p[i];  // update the line direction
			p[i]=pnew[i];        // update the starting point
		}

    #ifdef HACK_TO_DIAGNOSIS
      qnOut << its+1 << "   " << fret << endl;

      stringstream readIter;
      readIter << its+1 << "_";
      readIter << fret;
      funcd.print_sol(p, fret, outCif, string("Iter_")+readIter.str());
    #endif

    // find the coordinate for which |x_previous-x_new|/max(|x_new|,1) is biggest 
		test=0.0;
		for (Int i=0;i<n;i++) {
			temp=abs(xi[i])/MAX(abs(p[i]),1.0);
			if (temp > test) test=temp;
		}
    // return if this value is smaller than TOLX
		if (test < TOLX)
    {

      #ifdef HACK_TO_DIAGNOSIS
        qnOut << "\nClosed because movement in x small\n";
        qnOut.close();
        outCif.close();
        fileCounter++;
      #endif

			return;
    }

		for (Int i=0;i<n;i++) dg[i]=g[i];  // save the old gradient
		funcd.df(p,g);     // get a new gradient at point p

    // find the coordinate for which |g_new|*max(|x_new|,1) / max(f_new,1) is biggest
		test=0.0;
		den=MAX(fret,1.0);
		for (Int i=0;i<n;i++) {
			temp=abs(g[i])*MAX(abs(p[i]),1.0)/den;
			if (temp > test) test=temp;
		}
    // return if this value is smaller than gtol
		if (test < gtol)
    {
      #ifdef HACK_TO_DIAGNOSIS
        qnOut << "\nClosed because movement in x small\n";
        qnOut.close();
        outCif.close();
        fileCounter++;
      #endif

			return;
    }
		for (Int i=0;i<n;i++)
			dg[i]=g[i]-dg[i];
		for (Int i=0;i<n;i++) {
			hdg[i]=0.0;
			for (Int j=0;j<n;j++) hdg[i] += hessin[i][j]*dg[j];
		}
		fac=fae=sumdg=sumxi=0.0;
		for (Int i=0;i<n;i++) {
			fac += dg[i]*xi[i];
			fae += dg[i]*hdg[i];
			sumdg += SQR(dg[i]);
			sumxi += SQR(xi[i]);
		}
		if (fac > sqrt(EPS*sumdg*sumxi)) {
			fac=1.0/fac;
			fad=1.0/fae;
			for (Int i=0;i<n;i++) dg[i]=fac*xi[i]-fad*hdg[i];
			for (Int i=0;i<n;i++) {
				for (Int j=i;j<n;j++) {
					hessin[i][j] += fac*xi[i]*xi[j]
						-fad*hdg[i]*hdg[j]+fae*dg[i]*dg[j];
					hessin[j][i]=hessin[i][j];
				}
			}
		}
		for (Int i=0;i<n;i++) {
			xi[i]=0.0;
			for (Int j=0;j<n;j++) xi[i] -= hessin[i][j]*g[j];
		}
	}
	throw("too many iterations in dfpmin");
}

//template <class T>
struct Funcd 
{
	/*Doub EPS;
	T &func;
	Doub f;
	Funcd(T &funcc) : EPS(1.0e-8), func(funcc) {}*/
	virtual Doub operator() (VecDoub_I &x) = 0;
//	{
//		return f=func(x);
//	}

	virtual void df(VecDoub_I &x, VecDoub_O &df) = 0;
/*	{
		Int n=x.size();
		VecDoub xh=x;
		Doub fold=f;
		for (Int j=0;j<n;j++) {
			Doub temp=x[j];
			Doub h=EPS*abs(temp);
			if (h == 0.0) h=EPS;
			xh[j]=temp+h;
			h=xh[j]-temp;
			Doub fh=operator()(xh);
			xh[j]=temp;
			df[j]=(fh-fold)/h;
		}
	}*/
};

#endif 